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Optimal Multi-Robot Path Planning on Graphs: 
Complete Algorithms and Effective Heuristics 
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Abstract 

We study the problem of optimal multi-robot path planning on graphs (MPP) over four distinct minimization objectives: the 
makespan (last arrival time), the maximum (single-robot traveled) distance, the total arrival time, and the total distance. In a 
related paper and LaValle| ( |20l31 , we show that these objectives are distinct and NP-hard to optimize. In this work, we focus 
on efficiently algorithmic solutions for solving these optimal MPP problems. Toward this goal, we first establish a one-to-one 
solution mapping between MPP and network-flow. Based on this equivalence and integer linear programming (ILP), we design 
novel and complete algorithms for optimizing over each of the four objectives. In particular, our exact algorithm for computing 
optimal makespan solutions is a first such that is capable of solving extremely challenging problems with robot-vertex ratio as 
high as 100%. Then, we further improve the computational performance of these exact algorithms through the introduction of 
principled heuristics, at the expense of some optimality loss. The combination of ILP model based algorithms and the heuristics 
proves to be highly effective, allowing the computation of l.x-optimal solutions for problems containing hundreds of robots, 
densely populated in the environment, often in just seconds. 


I. Introduction 

In this paper, we study the problem of optimal multi-robot path planning on graphs (MPP), focusing on the design of 
complete algorithms and effective heuristics. In an MPP instance, the robots are uniquely labeled (i.e., distinguishable) and 
are confined to an arbitrary connected graph. The robots may move from a vertex to an adjacent one in one time step in the 
absence of collision, which may occur when two robots simultaneously move to the same vertex or along the same edge in 
different directions. A distinguishing feature of our MPP formulation is that we allow robots on fully occupied cycles to rotate 
synchronously. Such a formulation, more appropriate for multi-robot applications, has not been widely studied (except, e.g., 
|Standley| ( |2010| l; [Standley and Korf| ( |201![ )). Over the basic MPP formulation, we look at four commonly studied minimization 
objectives: the makespan (last arrival time), the maximum (single-robot traveled) distance, the total arrival time, and the 
total distance. We note that these global objectives have direct relevance toward real world multi-robot applications such as 


autonomous warehouse systems Wurman et al. (20081. For example, minimizing makespan is equivalent to minimizing the 


task completion time, whereas minimizing total distance is applicable to minimizing the fuel consumption of the entire fleet 
of robots. 

In a related paper Yu and LaValle (2015|l, we show that these objectives are all distinct and NP-hard to optimize, suggesting 


that efforts on solving optimal MPP should be directed at finding effective, near-optimal algorithms. In this work, we make an 
attempt toward this goal and propose a novel yet general framework for solving MPP optimally. By examining space and time 
dimensions jointly, we observe a one-to-one mapping between a solution for an MPP instance and that for a multi-commodity 
flow problem derived from the MPP problem. Based on the equivalence relationship, we can then translate the MPP problem 
into an integer linear programming (ILP) model, which can be subsequently solved using an ILP solver. The generality of 
ILP allows the encoding of all four objectives to yield complete algorithms for optimizing these objectives. From here, we 
take a further step and introduce several heuristics to boost the algorithmic performance at a slight loss of solution optimality. 
Our method is especially effective in computing near-optimal minimum makespan solutions, capable of computing l.x-optimal 
solutions for hundreds of robots, densely populated on the underlying graph, often in just seconds. 


Related work. Multi-robot path planning problems, in its many formulations, have been actively studied for decades |Erdmann 

and Lozano-Perez ( 

1986)1; 

Zelinsky 

( 1992)1 ; LaValle and Hutchinson 

1998)1; Silver)()2005)l; Kloder and Hutchinson 

()2006 1; Ryan 

(|2008i; Jansen and Sturtevant (2008 1 ; Surynek (2009 1 ; Luna and Bekris (201 1)1; |Standley and Korf|(|201 l|l; van 

den Berg and 

Overmars (2005)1; 

van den Berg et al. (2009)1; Solovey and Halperin (2012)1; )Yu and LaValle (2013a)l; Turpin et al. (2014)1; 

Solovey et al. P0151l. As 

a universal subroutine, collision-free path planning for multiple robots finds applications in tasks 


spanning assembly Halperin et al. 

(2000]l; iNii^ ( 

1992)1, evacuation Rodriguez and Amato 

()2010)l, formation control jBalch 

and Arkin 

(| 1 998)1 ; IPoduri and Suk] 

latme (]2004)l; ]S 

tucker et al. ( 

20o7)i; jSmith et al.) ()2000)l; 

Tanner et al.)()2004|, localization 

Fox et al. 

()2000)l, micro droplet manipulation Ding et al.) ()200 

)l; Griffith and Akella (2005)1, object transportation Mataric 

et al.)(1995i; Rus et al. ()1995l, search and rescue Jennings et al. (1997)1, and so on. We refer the readers to Choset et al. 


le (20061 and the references therein for a more comprehensive review on the general subject 


of multi-robot path and motion planning. 
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The algorithmic study of graph-based multi-robot path planning problems, the focus of this work, can be traced to as least 
1879 Story ( 1879| l, in which Story makes the observation that the feasibility of the 15-puzzle |Loyd| ( [1959[ ) is decided by the 
parity of the game. The 15-puzzle is a restricted MPP instance moving 15 labeled game pieces on a 4 x 4 grid, from some 
initial configuration to some goal configuration. The restriction is that only a single game piece near the only empty vertex 
may move to the empty vertex in a step-multiple mobile robots could move simultaneously. A generalization of the 15-puzzle 
is introduced in Wilson ( |1974 1 , extending the problem from 15 game pieces on a 4 x 4 grid to n — 1 labeled pebbles on an 
n-vertex, 2-connected graph. It is shown, together with an implied algorithm, that an instance is always feasible if the graph 
is non-bipartite. When the graph is bipartite (such as the 15-puzzle), all pebble configurations are split into two groups of 
equal size such that any two configurations in the same group form a feasible instance. A further generalization is introduced 
(1984|, allowing p < n pebbles on a graph with n vertices. For this problem, an 0{n^) algorithm is 


in 


Kornhauser et al. 


provided to solve an instance or decide that the instance is infeasible. 

As computer games and multi-robot systems gain popularity, concurrent movements are introduced and pebbles are replaced 
with robots (or agents). On the feasibility side, the MPP problem studied in this paper is shown to be solvable also in 0{n^) 
time Yu and Rus ( 2014| l. To distinguish the formulations, we denote the formulation that does not allow cyclic rotations of 
robots along fully occupied cycles as cycle-free MPP. Until recently, the majority of algorithmic study on MPP is on the cycle- 
free case. Since the problem is shown to be tractable Kornhauser et al. ( 1984| l, most algorithmic study of cycle-free MPP put 
some emphasis on optimality. Through the clever use of primitive operations, algorithms from [Ryan] ( |2008| l; |Surynek| ( |2009| l; 
|Luna and Bekris|p011| l; [Sajid et ah] ( |2012| l; |de Wilde et al.| ( |2014| l could quickly solve difficult problems with some form of 
completeness guarantees. These algorithms do not have optimality guarantees but the produced solutions are often of much 

(1984i. For more discussion and references on sub-optimal 


better quality than the 0{rt') bound given by 

Kornhauser et al. 

00 

o^ 

methods, see Sajid et al.l (|2012|); de Wilde et a 

1. ( 

201411. 


Pushing more toward the optimality side, most algorithmic results explore ways to limit the exponential search space growth 
induced by multiple robots. One of the first such algorithms. Local Repair A* (LRA*) |Zelinsky| ( |1992| l, plans robot paths 
simultaneously and performs local repairs when conflicts arise. Focusing on fixing the (locality) shortcomings of LRA*, 
Windowed Hierarchical Cooperative A* (WHCA*) proposes using a space-time window to allow more choices for resolving 
local conflicts while simultaneously limiting the search space size |Silv^ (2005 1 . Iterative deepening technique is shown to be 
applicable to MPP problems in Sharon et al. ( 2012| l. A technique called sub-dimensional expansion is shown to perform well 
in complex environment Wagner and Choset ( 201 1[); th e robot density is however relatively low (104 cells per robot according 
to the paper). In Standley ( 2010| l; Standley and Korf ( 201 l) l, instead of agnostic dissection of an instance, the natural idea 
of independence detection (ID) is explored to only consider multiple robots jointly (the source of exponential search space 
growth) as necessary. With operator decomposition (OD) that treats each legal move as an “operator”, the authors produced 
algorithms (ID,ODh-ID, and related variants) that prove to be quite effective in computing total time- or distance-optimal 
solutions. We point out that ID and ODh-ID have support for handling cycles (i.e., they apply to MPP instead of cycle-free 
MPP). Algorithms designed for minimizing makespan have also been attempted, e.g., |Surynek| ( |201^ , but the solution quality 
degrades rapidly as the robot-vertex ratio increases. 

Many approaches have also been proposed for solving multi-robot path planning problems in the continuous domain. A 


representative method called velocity-obstacles Kant and Zucker ( 

C) 

00 

o^ 

v^an den Berg et al. (12008 201111 explicitly examines 

velocity-time space for coordinating robot motions. In Griffith and Akell 

a (2005|l, mixed integer programming (MIP) models 


are employed to encode the robot interactions. A method based on the space-time perspective, similar to ours, is explored in 
1. Karamouzas (20121. In|Peasgood et al.' (2008|l, an A*-based search is performed over a discrete roadmap abstracted from 


the continuous environment. In Solovey et al. (20141, discrete-RRT (d-RRT) is proposed for the efficient search of multi-robot 


roadmaps. Algorithms for discrete MPP, cycle-free or not, have also helped solving continuous problems Solovey and Halperin 
P012| l; IKrontiris et al.| ( |2012| l. 

Contributions. We study the optimal MPP formulation allowing up to n robots on a n-vertex connected graph, which 
we believe is better suited for multi-robot applications. The formulation is not widely studied, perhaps due to the inherent 
difficulty in handling cyclic rotations of robots. Besides the novelty of the problem, this work brings several algorithmic 
contributions. First, based on the equivalence relationship between MPP and network flow, we establish a general and novel 
solution framework allowing the compact encoding of optimal MPP problems using ILP models. We show that the framework 
readily produces complete algorithms for minimizing the makespan (last arrival time), the maximum (single-robot traveled) 
distance, the total arrival time, and the total distance, which are perhaps four most common global objectives for MPP. The 
resulting algorithms, in particular the one for computing minimum makespan, are highly effective in computing challenging 
problems instances with robot-vertex ratio up to 100%. Second, we introduce several principled heuristics, in particular a k-way 
split heuristic that divides an MPP instance over the time horizon, to give the exact algorithms a sizable performance boost at 
the expense of some solution optimality loss. With these heuristics, we are able to extend our algorithms to tackle problems 
with several hundred robots that are extremely densely populated, while at the same time maintain 1.x solution optimality. Last 
but not least, our successful exploit of ILP to attack optimal MPP shows that integer linear programming method is competitive 
with direct search methods in this problem domain, especially when the number of robots becomes large. This is surprising 
because ILP solvers are not designed specifically for MPP. 
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The rest of the paper is organized as follows. In Section |II] we define the optimal MPP problems and provide background 


material on network flow. We then establish the equivalence relationship between MPP and multi-commodity flow in Section III 


We derive complete algorithms in Section IV based on the equivalence and continue to describe the performance boosting 
heuristics in Section [V] We evaluate the algorithms in Section |VI| and conclude in Section VII This paper is partly based 


Yu and LaValle 


( |2013a|b|c| )n In comparison 


to 


Yu and LaValle 


(2013a|b I, besides demonstrating significantly improved 


computational performance due to the addition of the fc-way split heuristic, we have substantially extended the generality of 
our ILP-based algorithmic framework, which now supports all common global time- and distance-based objectives. 


II. Preliminaries 

We now define MPP and the optimality objectives studied in this paper. Following the problem statements, we gave a brief 
review on network flow. 


A. Multi-Robot Path Planning on Graphs 

Let G = (y, E) be a connected, undirected, simple graph, with V = {vi} being the vertex set and E = {(ui, vj)} the edge 
set. Let R = {ri,..., r„} be a set of n robots. The robots move at discrete time steps (i.e., at f = 0,1,...). At time step t = 0, 
each robot occupies a distinct vertex of G. In general, at any time step f = 0,1,..., the robots assume a conflguration that is 
an injective map from R to V. The start {initial) and goal configurations of the robots are denoted as xi and xq, respectively. 
Fig. [2 a) shows a possible conflguration of 9 robots on a 3 x 3 grid graph. Fig. [^b) shows a possible goal conflguration, in 
which the robots are ordered based on row-major ordering 


(a) (b) 

Fig. 1. a) A 9-puzzle problem, b) The desired goal configuration. 




During a discrete time step, each robot may either remain stationary or move to an adjacent vertex. To formally describe 
a plan, let a scheduled path be a map pi : —)• V, in which := N U {0}. A scheduled path pi is feasible if it satisfies 

the following properties: 1) Pi(0) = xiirf). 2) For each i, there exists a smallest L € such that pflti) = xc{ri). 3) For 
any t > U, pjf) = xairi). 4) For any 0 < f < U, {pflt),pi{t -\- 1)) e E or pflt) = pflt + 1) (if Pi{f) = pflt -f 1), robot r* 
stays at vertex pflt) between the time steps t and t -f 1). We say that two paths Pi,Pj are in collision if there exists k G Z+ 
such that Pi{t) = Pj(t) (meet collision) or {pi{t),pi{t -\- 1)) = {pj{t -f l),pj{t)) (head-on collision). As an illustration. Fig. 
shows the feasible and infeasible moves for two robots during a single time step The multi-robot path planning on graph 
(MPP) problem is defined as follows. 

(a) (b) (c) 

Fig. 2. Some feasible and infeasible moves for two robots, a) A feasible synchronous move, b) An infeasible synchronous move in which two robot collide 
“head-on”, c) An infeasible synchronous move in which two robots “meet” at a vertex. 


Problem 1 (Multi-robot Path Planning on Graphs) Given a A-tuple {G, R, xi, xq), find a set of paths P = {pi,... ,p„} 
such that Pi's are feasible paths for respective robots Vi ’s and no two paths pi , pj are in collision. 


For example. Fig. fTta) and Fig. fTtb) define an MPP problem on the 3x3 grid. We call this particular problem the 9-puzzle 

de| ' ' ' 


problem (a variant of the 15-puzzle |Ratner and Warmuth (1990i), which readily generalizes to W^-puzzles. 

Remark. With a few exceptions {e.g., Standley and Korf ( 201 l| l), most existing studies on discrete multi-robot path planning 
problems require empty vertices as swap spaces. In these formulations, in a time step, a non-intersecting chain of robots 
may move simultaneously only if the head of the chain is moving into a previously unoccupied vertex. In contrast, our MPP 


Yu and LaValle 


1 2013c 


is a preliminary poster presentation. 


^In this paper, we generally use shaded discs to mark start locations of robots and discs without shading for goal locations. 

^We assume that the graph G only allows “meet” or “head-on” collisions. The assumption is mild. For example, a (arbitrary dimensional) grid with unit 
edge distance is such a graph for robots of with radii of no more than v^/4. 
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formulation allows synchronized rotations of robots along fully occupied cycles (see, e.g.. Fig. [T]and[^. This implies that even 
when the number of robots equals the number of vertices, robots can still move on disjoint cycles. We note that MPP can be 
solved in polynomial time with feasibility test taking only linear time Yu and Rus (|2014|l. A 


B. Optimal Formulations 

Let P = {pi,... ,p„} be an arbitrary feasible solution to some fixed MPP instance. For a path pi G P, let len{pi) denote the 
length of the path pi, which is increased by one each time when the robot passes an edge. A robot, following a path pi, may 
visit the same edge multiple times. Recall that L denotes the arrival time of robot In the study of optimal MPP formulations, 
we examine four common objectives with two focusing on time optimality and two focusing on distance optimality. Below, 
each objective is formally defined. 

Objective 1 (Makespan) Compute a path set P that minimizes maxi<i<„ ti. 

Objective 2 (Maximum Distance) Compute a path set P that minimizes maxi<i<„ len{pi). 

Objective 3 (Total Arrival Time) Compute a path set P that minimizes 

Objective 4 (Total Distance) Compute a path set P that minimizes X]r=i l^ti{pi). 

The intuitive meaning of these objectives is clear. Fig. [^illustrates the four-step minimum makespan solution to the 9-puzzle 
problem from Fig. [T] The solution is optimal as it takes at least four steps for robot 9 to move to its goal. 



Fig. 3. 


A 4-step solution from our algorithm. The directed edges show the moving directions of the robots at the tail of the edges. 


C. Network Flow Review 

A network JV = {G,ci,C 2 , S) consists of a directed graph G = {V,E) with ci,C 2 : E —>■ as the maps specifying 

capacities and costs over edges, respectively, and S' C V as the set of sources and sinks. Let S = S’*" U S~, with S^ denoting 
the set of source vertices, S~ denoting the set of sink vertices, and S’*" n S~ = 0. For a vertex v G V, let (5+(u) (resp., 
(5“(t;)) denote the set of edges of G going to (resp., leaving) v. A feasible (static) S^, S”-flow on this network Af is a map 
f : E —>■ that satisfies edge capacity constraints, 

Ve e E, /(e) < ci(e), (1) 

the flow conservation constraints at non terminal vertices, 

WvGV\S, Y, fie) - Y /(e) =0, (2) 

eG(5+(D) 

and the flow conservation constraints at terminal vertices, 

Pif) = E ( E /(e) - E /(e)) 

ves+ ee6-{v) ees+{v) .... 

= E ( E - E «'))■ 

v^S~ eG(5+(t') eG(5“(i') 

The quantity E{f) is called the value of the flow /. The classic (single-commodity) maximum flow problem asks the following 
question: Given a network N, what is the maximum F{f) that can be pushed through the network? The minimum cost maximum 
flow problem further requires the flow to have a minimum total cost among all maximum flows. That is, we want to find a 
flow among all maximum flows that also minimizes the quantity 

E'= 2 (e)-/(e). (4) 

eeE 

The network flow formulation described so far only considers a single commodity, corresponding to all robots being 
interchangeable. For general MPP formulations, the robots are distinct and must be treated as different commodities. Such 
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problems can be captured with the Multi-commodity flow or simply multiflow. Instead of having a single flow function /, we 
have a flow function fi for each commodity i. The constraints and become 

^ /j(e) < ci{e), (5) 

i 

Y. - E Me) = 0, (6) 

eG(5''"(D) eG<5~(T') 


Vi, 


E ( E f.M 

'uGS+ 6^5 ("u) 

E ( E 


vGS eG<5+(t)) 


- E 

eG(5+ (v) 

- E /*("))• 

eGS~ ( tj ) 


(7) 


Maximum flow and minimum cost flow problems may also be posed under a multiflow setup; we omit the details. Our review 
on network flow only touches aspects pertinent to the current work; for a thorough coverage on the subject of network flow, 
Aronson ( 1989|l; Ahuja et al. ( 1993|l and the references therein. Note that the multiflow model stated here is sometimes 


see 


also referred to as integer multiflow because fi must have integer values. 


III. From Multi-Robot Path Planning to Multiflow 

A close algorithmic connection exists between optimal MPP and network flow problems. Maximum (single commodity) flow 
problems generally admit efficient (low-degree polynomial time) algorithmic solutions Ahuja et al. ( 1993| l, whereas maximum 
multiflow is a well-known NP-hard problem, difficult to even approximate Andrews and Zhang ( 2005| l. Mirroring the disparity 
between single- and multi-commodity flows, in the domain of MPP problems, if there is a single group of interchangeable 
robots (here, for a group, it does not matter which robot goes to which goal as long as all goal locations assigned to the 
group are occupied by robots from the same group), then many optimal formulations admit polynomial time algorithms |Yu and] 
LaValle (2013aI. However, as soon as a single group of robots splits into two or more groups, finding optimal paths for these 
robots become intractable |Yu and LaValle ( |2015| l. The apparent similarity between optimal MPP and multiflow is perhaps best 
explained through a graph-based reduction from MPP problems to network flow problems. The reduction will also form the 
basis of our algorithmic solution. 

To describe the reduction, we use as an example the undirected graph G in Fig. Qa), with start vertices {sj *"},i = 1,2 and 
goal vertices = 1,2. An instance of MPP is given by (G, {ri, r2},x/ : i—>■ sf,xo : i—>■ s~). We will reduce the 

problem to a network flow problem JV = {G',ci,C 2 , S). The reduction proceeds by constructing a network that is a time- 


u{t)' u{t-\-l) 

\ / 

I 

/ 

• • 

v{t)' u(f+l) 

(b) 

Fig. 4. a) A simple G. b) A merge-split gadget for splitting an undirected edge through time steps, for enforcing the head-on collision constraint. 

expanded version of the graph G, which then allows the explicit consideration of the interactions among the robots over space 
and time. To carry out this expansion, a time horizon must first be decided. For different optimality objectives, the expansion 
time horizon, some natural number T, is generally different; for now we assume that T is fixed. 

To begin building the network, we create 2T -f 1 copies of G’s vertices, with indices 0,1,1',..., as shown in Fig. For 
each vertex v G G, denote these copies u(0) = t;(0)', w(l), t;(l)', u(2),..., u(r)'. For each edge {u,v) £ G and time steps 
t,t-\-l, 0 < t < T, the merge-split gadget shown in Fig. |^b) is added between u{t)', v{ty and u{t -f 1), v{t -f 1) (arrows from 
the gadget are omitted from Fig. since they are small). For the gadget, we assign unit capacity to all edges, unit cost to the 
horizontal middle edge, and zero cost to the other four edges. The merge-split gadget ensures that two robots cannot travel 
in opposite directions on an edge in the same time step, which prevents head-on collision between two robots. To finish the 
construction of Fig. for each vertex v G G, we add one edge between every two successive copies (i.e., we add the edges 
(t;(0), u(l)), (ujl), t;(l)'),..., (ujT), u(r)')). These correspond to the green and blue edges in Fig. For all green edges, 
we assign them unit capacity and cost; for all blue edges, we assign them unit capacity and zero cost. The green edges allow 



’2 4 


(a) 
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Fig. 5. The time-expanded network with an expansion time horizon of T = 2 over the base graph Fig. Ba). 


robots to stay at a vertex during a time step whereas blue edges ensure that each vertex holds at most one robot, enforcing 
the meet collision constraint. 

Fig. [^(with the exception of edges e\ and 62 , which will become relevant shortly), the time-expanded network, is the desired 
G. For the set S' = S'"*" U S“, we may simply let S+ = {t^(0) : v G {s^}} and S~ = {u(T)' : v G {s“}}. Af = {G, ci, C 2 , S) 
is now complete; we have reduced MPP to an integer multiflow problem on JV with each robot from i? as a single type of 
commodity. 

Theorem 1 Let (G, R, xj, xa) be an MPP instance. There is a bijection between its solution set (with a maximum number of 
time steps up to T) and the integer maximum multiflow solutions of flow value n on the time-expanded network Af constructed 
from (G,R,Xi^xg) with T time steps. 

Proof. Injectivity. Assume that P = {pi,... ,p„} (n is the number of robots) is a solution to an instance of MPP. For each 
Pi and every time step t = 0,... ,T, we mark the copy of pflf) and pflt)' (recall that pfli) corresponds to a vertex of G) 
at time step t in the time-expanded graph G. Connecting these vertices of G sequentially (there is a unique way to do this) 
yields one unit of flow fi on Af (after connecting to appropriate source and sink vertices in which is trivial). It is 

straightforward to see that if two paths Pi^Pj are not in collision, then the corresponding flows fi, fj on Af are vertex disjoint 
paths and therefore do not violate any flow constraint. Since any two paths in P are not in collision, the corresponding set of 
flows {fi,..., fn} is feasible and maximal on Af. 

Surjectivity. Assume that {/i,..., /„} is an integer maximum multiflow on the network Af with \ fi\ = 1 for all i’s. First 
we establish that any pair of flows fj are vertex disjoint. To see this, we note that fj (both are unit flows) cannot share 
the same source or sink vertices due to the unit capacity structure of Af enforced by the blue edges. If fi,fj share some 
non-sink vertex v at time step t > 0, both flows then must pass through the same blue edge (see Fig. |4]jr)) with v being either 
the head or tail vertex, which is not possible. Thus, fj are vertex disjoint on Af. We can readily convert each flow to a 
corresponding path pi (after deleting extra source vertex, sink vertices, vertices in the middle of the gadgets, and tail vertices 
of blue edges) with the guarantee that no Pi,Pj will collide due to a meet collision. By construction of Af, the gadget we used 
ensures that a head-on collision is also impossible. The set {pi,... ,pn} is then a solution to MPP. □ 


Remark. A multiflow problem can be reduced to an MPP problem (on a directed graph) as well. In particular, if all edges in 
the network have unit capacity and there is a single unit of each type of commodity, then a multiflow problem is a multi-robot 
path planning problem, often known as the edge disjoint path problem Robertson and Seymour (1995i. A 


IV. Complete, Integer Linear Programming-Based Algorithms eor Optimal MPP Problems 

Because optimizing MPP solutions over Objectives [T]|2 are computationally intractable, reducing MPP to multiflow problems 
does not make these optimal MPP problems any easier. However, with a network flow formulation (see Secion [lL^ ), it becomes 
possible to establish integer linear programming (ILP) models for optimal MPP formulations. These ILP models can then be 
solved with powerful linear programming packages. In comparison to A*-based algorithms augmented with heuristics which 
often target important but a limit set of problem structures, ILP-based algorithms proposed here are agnostic to specific problem 
structures. As such, ILP-based algorithms appear more capable at attacking a wider range of MPP problems and in particular 
difficult MPP instances in which the robot-vertex ratio is high. In this section, we build ILP models for each of Objectives [T]|^ 

A. Minimizing the Makespan 

A minimum makespan solution to an MPP instance / = {G, R,Xi,xg) can be computed using a maximum multiflow 
formulation. Fixing a time span T, let N = {G, Ci, C 2 , S) be the time-expanded network for I, a set of n loopback edges are 
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added to G' by connecting each pair of corresponding start and goal vertices in S, from the goal to the start. We use e^ ’s 
to denote edges of G' and let the n loopback edges take the hrst n indices, with , 1 < j < n being the edge connecting 
the goal vertex of rj to the start vertex of Vj. For example, for the G' in Fig. ei and 62 are the loopback edges for ri 
and r 2 , respectively. The lookback edges have unit capacities and zero costs. Next, for each edge ej G G' (including the 
lookback edges), create n binary variables cci j ,..., Xn,j corresponding to the flow through that edge, one for each robot. The 
is, Xij = 1 if and only if robot passes through Cj in G'. The variables Xij's must satisfy two edge capacity constraints and 
one flow conservation constraint, 

n 

i—1 

yi<ij <n, j, Xij = 0, 


The objective function is 


Wv G G' and 1 < i < n, E XiJ — E Xij. 

ejGS+{v) ejGS-{v) 

max Xi^i. 

l<i<n 


(9) 

( 10 ) 


For each fixed T, a solution to the above ILP model with an objective equaling n means that a feasible solution to MPP is 
found. We are to hnd the minimal T that yields such a feasible solution. To do this, start with T being the maximum over all 
robots the shortest possible path length for each robot, ignoring all other robots. An ILP model for this T is then built and 
tested for a feasible solution. If the model is not feasible, T is increased and the procedure repeated. The hrst feasible T is the 
optimal T. Once a feasible model with the minimum expansion time horizon T is found, the robots’ paths can be extracted 
based on the proof of TheoremThe algorithm is in fact complete: Since the problem is discrete, there is only a hnite number 
of possible states. Therefore, for some suffici ently large T , there must either be a feasible solution or we can pronounce that 
none can exist (T = 0(|Fp) is big enough 


Yu and Rus 


(2014 1 ). Denoting the resulting algorithm as MinMakespan 


have shown the following. 


Proposition 2 Algorithm MinMakeSPAN is a complete algorithm for finding minimum makespan solutions for MPP. 


B. Minimizing the Maximum Single-Robot Traveled Distance 

For minimizing the maximum distance traveled by any single robot, the network and variable creation remains the same as 
the minimum makespan setup; constraints ([^ and (|^ remain unchanged. Because we want to send all robots to their goal, a 
maximum flow may be forced through the additional constraint 

VI < 7 < n, Xi^i = 1. (11) 

To encode the min-max objective function, we introduce an additional integer variable Xmax and add the constraint 

^ ^ C2(Cj) * Xi j ^ XjYiax (12) 

ej^G' ,j>n 

for all 1 < 7 < n. For a fixed i, the left side of ( [T^ represents the distance traveled by robot r^. The objective function is 
then simply 

min Xmax- (13) 

Denoting the algorithm as MinMaxDist, we have 

Proposition 3 Algorithm MinMaxDist is a complete algorithm for finding solution to MPP that minimize the maximum 
distance traveled by a single robot. 

Remark. Assuming that we have a minimum makespan solution, we can better bound the time horizon T needed for 
MinMaxDist. Let tmm be the minimum makespan computed by MinMakespan, setting T = ntmin is then sufficient 
for finding the min Xmax- To see that this is true, we first note that min Xmax < tmin because the minimum makespan solution 
requires robots to to synchronize their moves, which may force some robots to travel unnecessarily. Then, the n robots cannot 
move a total distance of more than n minx max < ntmin because no robot may travel more than minx max edges and therefore 
cannot use more than n min Xmax time in total. A 





C. Minimizing the Total Arrival Time 


The ILP model for minimizing the total arrival time is more involved in the way the objective function is constructed. First, the 
network and all variables from the minimum makespan ILP-model are inherited. We also inherit constraints ([^, (0, and ( [TT] l. 
To represent the objective function, for each time step \ <t <T and each v := xc{ri), 1 < i < n, we create a binary variable 
y*. Then, we give new indicies to some variables that are already created. Recall that for each edge ej = {v{t),v{ty) G G' 
(e.g., the four extra bold blue edges in Fig. 1^, a variable Xij is created. There are a total of nT such variables. Here, we give 
these variable a second index xj. That is, cq is the binary variable indicating whether edge {v{t),v(ty) G G', v := xairi), is 
used by robot r^. 

Given a network with a fixed T, if constraints and o can be satisfied, then there is a feasible solution to the 

original MPP problem. In this case, xf = 1 for 1 < i < n. We let yf = xj. Then, each y\, 1 < f < T is defined recursively 
over x\ and as 


vl > J/i yl < y^ 


y\ < x\. 


(14) 


The constraint ( [T4l i effectively performs the logical and operation over x\ and yl'^^ and stores the result in y\. In the end, 
the smallest t for which y* = 1 is the time robot reaches its goal (and stops). Therefore, for each 1 < i < n, ELi yl is the 
number of time steps from the time arrives at its goal until time T. T — Y^=i yl is then the time spent by r^. To minimize 
the total arrival time, the objective function can be expressed as 


nT- yl 

l<i<n,l<t<T 


(15) 


Denoting the resulting algorithm (with a sufficiently large T) as MinTotalTime, we have 


Proposition 4 Algorithm MinTotalTime is a complete algorithm for finding minimum total arrival time solutions for MPP. 


D. Minimizing the Total Distance 

From the ILP model for minimizing the maximum distance, we only need to change the objective function for computing 
a minimum total distance solution. We do not need the variable Xmax and simply update the objective function to 

min E C2{ej)-Xij. (16) 

ej^G' ,j>n, l<2<n 

Denoting the algorithm as MinTotalDist, we have 

Proposition 5 Algorithm MinTotalDist is a complete algorithm for finding minimum total distance solutions for MPP. 
Again, T = ntmin is sufficient for building a network that contains a minimum total distance solution, if one exists. 


V. Heuristics eor Eeeective Computation of Near-Optimal Solutions 

In Section IV in constructing the ILP models for Objectives T]|4 our goal is to show the universal applicability of the 
network flow model {e.g.. Fig. toward many different optimal MPP formulations. The resulting ILP-based algorithms are 
complete and always produce the optimal solution in principle. As will be shown in Section [V^ such algorithms perform very 
well in handling relatively small but extremely challenging problems. Nevertheless, as the problem size grows (i.e., as the 
graph G and the number of robots n get larger), the computation time needed by ILP solvers grows rapidly. From a practical 
point of view, it may be far more desirable to quickly compute a good quality but sub-optimal solution than waiting forever 
for the optimal solution. In this section, we introduce several heuristics to accomplish just that, with a particular focus on 
computing solutions with minimum makespan. 


A. Building More Compact Models 

To get the best performance out of a solver, it is beneficial to have a lean model (i.e., fewer columns and rows). So far, our 
focus has been to provide a general network flow based framework so that the ILP models can be easily built. When it comes 
to translating the models to an ILP solver, they can be further simplified. The heuristics discussed in this subsection aim at 
making the representation of the constraints ([^ and (|^ more compact in the resulting ILP models. As such, they apply to all 
the optimality objectives. 
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• Si 


Fig. 6. A more compact representation of the network flow graph from Fig. 


u{t) u{t-\-l) 



Fig. 7. The simplified merge-split gadget for enforcing the head-on collision constraint. 


Better encoding of the collision constraints: Recall that in building the network flow model (e.g.. Fig. |^, we used a merge- 
split gadget (Fig. Qb)) for enforcing the head-on collision constraint and extra time steps {e.g., the blue edges in Fig. for 
avoiding meet collisions. When we translate this into linear constraints, these structures can be simplified to yield the more 
compact structure illustrated in Fig. 

In the newer structure, each merge-split gadget now has two edges instead of five. Also, the blue edges are removed. The 
updated gadget for an edge (u, v) G E between time steps t and f -f 1 is shown in Fig. (note that due to the removal of 
the blue edges, vertices such as v{ty is no longer needed). That is, instead of five, only two variables are needed for each 
robot. Denoting these binary variables as Xi^(^u{t),v(t+i)) for a robot r^, the head-on collision constraint for 

a single gadget can be readily encoded as 

n n 

Xi,{v(t),u{t-i-l)) < 1- (17) 

i=l i=l 

Then, to enforce the meet constraint, for example at a vertex v{t), we simply require that at most one outgoing edge from 
v{t) may be used, i.e., 

Y. ^^7 < 1- (18) 

ej {v{t)),l<i<n 

Overall, the newer ILP model is roughly half of the size of the original model. 

Reachability analysis: In the time-expanded graph, there are redundant binary (edge) variables that can never be true because 
some edges are never reachable. For example, in Fig. at f = 0, the only outgoing edges that can possible be used are those 
originates from and sj. The rest can be safely removed. In general, for each robot r^, based on its reachability from its 
start vertex and to its goal vertex, a sizable number of binary variables Xij's can be deleted. 

B. Divide-and-Conquer Over Time Domain 

In evaluating the ILP model-based algorithm for optimal makespan computation, we observe that the ILP solver running 
time appears to grow exponentially as the size of the model grows. This prevents the algorithm from performing well over 
instances with more than a few tens of robots. The observation, while hampering the effectiveness of the exact algorithm, 
turns out to offer a useful insight toward a highly effective heuristic. We find that, when the overall size of the ILP model 
is relatively small and the robot-vertex ratio is not approaching 1, even when there are a large number of robots, the model 
usually does not present much challenge for ILP solvers. To apply the ILP model-based method to more challenging problems 
(e.g., solving problems with hundreds of robots quickly), we simply limit the size of the individual ILP model fed to the 
solver. One way to achieve this is through divide-and-conquer over the time domain. We use a simple example (see Fig. to 
illustrate the idea. 
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Fig. 8. a) A simple two-robot problem, b) The time-divided instances. 


In Fig. I^a), we have a simple planning problem for two robots on a 3 x 3 grid. To execute the heuristic, we first compute 
a shortest path for every pair of start and goal locations. In this case, we get the orange and green paths for robots 1 and 
2, respectively. Then, if we decide to split the problem into two smaller problems, for each of the paths, it is split into two 
(generally) equal length pieces and the middle node is set as the intermediate goal. In our example, we may do this for robot 
1 by setting the intermediate goal location at (1,1) from the top-left corner (the brown disc labeled 1 in Fig. |^b)). For robot 
2, because the middle location coincides with that of robot 1, we pick an alternative location that is not already occupied as 
the intermediate goal for robot 2, in this case (2, 2) from the top-left corner. The intermediate goals for the first instance will 
also serve as the start locations of the second instance. This yields two child instances with both requiring a time expansion 
with 2 steps each, effectively making the individual ILP model roughly half the size of the original one, which requires a time 
expansion with 4 steps. In general, we may divide a problem into arbitrarily many smaller instances in the time domain. 

If a problem is divided in this manner to k sub problems, we call the resulting heuristic a k-way split. Because the 
division is over time, there is no interaction between the individual, smaller instances. Once we obtain the solution for each 
child instance, the solutions can be glued together by concatenating the results. In practice, this simple heuristic dramatically 
improves algorithm performance without heavy negative impact on path optimality in terms of makespan; we observe a 
consistent speedup in computational experiments. 

Remark. The k-way split heuristic, by design, is particularly suitable for the makespan objective. This is the case due to the 
additive nature of the makespan objective over the split sub-problems. Besides makespan, the heuristic also applies to distance 
objectives {i.e.. Objectives andquite well, as long as the time horizon required for finding distance optimal solution does 
not differ greatly from the time horizon required for minimum makespan solution. The heuristic does not directly apply to 
Objective [^because total time is not additive over the split sub-problems. As an example, suppose that a 2-way split is carried 
out with each sub-problem having a time horizon of T /2. If a robot does not move in the solution to the first sub-problem 
{i.e., 0 <t < T/2), it contributes 0 to the total distance. However, if r; moves even a single step in the solution to the second 
sub-problem (i.e., T/2 <t < T), then will contribute at least T/2 to the total arrival time. Nevertheless, k-way split is still 
helpful in this case as we may use it to quickly compute an initial T for performing the time expansion. A 


VI. Experimental Evaluation 


In this section, we evaluate the performance of our optimal and near-optimal MPP algorithms with an emphasis on 
MinMakespan. Our performance evaluation covers a broad spectrum of typical problem settings. Eor each setting, we 
push the limit on the robot-vertex ratio-to as high as 100%. To the best of our knowledge, the majority of the settings with 
high robot-vertex ratio have never been attempted with much success prior to our study. 

When applicable, we also compare our results with the state-of-the-art found in the literature. In particular, we have examined 
ODh-ID, id, WHCA*, and COBOPT Surynek| ( 2012| l, among others. ODh-ID and ID support cycles whereas WHCA*appears 
to be designed for cycle-free MPP as it could not solve any 9-puzzle. The problem definition for C OB OPT suggests it solves MPP 
but it employs a cycle-free subroutine for finding feasible solutions. ODh-ID, ID, and WHCA* are designed for optimizing 
total time and total distance optimal objectives, and do not naturally extend to makespan computation. However, the associated 
makespan produced by these algorithms are usually of good quality. Among these three, our experiments show that ID-based 
anytime algorithm is the most versatile due to its IDA*-like incremental structure. On the other hand, ODh-ID and WHCA* 
do not scale well when the robot-vertex ratio goes beyond 10-20%, depending on the particular problem setting]^ COBOPT is 
designed for makespan computation. 

We implemented all algorithms (MinMakespan, MinMaxDist, MinTotalTime, and MinTotalDist) in the Java 
programming language. We take advantage of multi-core CPUs when the k-way split heuristic is being used. Also, Gurobi 
Gurobi Optimization ( 2014| l, the ILP solver used in our implementation, can engage multiple cores automatically for hard 
problems. We ran all our tests on a MacBook Pro laptop computer (Intel Core i7-4850HQ, 16GB memory). We thank Trevor 
Standley for sharing the C code implementing ODh-ID, ID, and WHCA*, among others. We modified (the original code 


^Some of these algorithms were evaluated without requiring all robots reach their goals. For example, in the WHCA* work |Silver| jlOOs) , if an instance 
with n robots is solved for p < n robots, the problem is counted as partially solved. We require each instance to be fully solved to be counted as a success. 
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supports only 32 x 32 grid) and compiled the code as a 64-bit windows executable under MSVC 2010 with all speed optimization 
flags turned on. The comparison to COBOPT uses the result provided in Surynek P012| l, which covers only 8x8 and 16 x 16 
grids. 


A. Performance of MinMakespan and k-way Split Heuristic 

We begin our experimental evaluation focusing on the MinMakespan algorithm and the k-way split heuristic. For this 
purpose, we use as the based graph a 24 x 18 grid with varying number of vertices (0-25%) removed to simulate obstacles. 
The connectivity of the graph is always maintained. We note that with 25% vertices removed, the graph is already sparsely 
connected at some places (see e.g., Fig. [^, making solving these problems optimally a challenging task. In presenting the 



Fig. 9. A typical 24 x 18 grid instance with 25% vertices removed to simulate obstacles and 50 start and goal locations. Note that the connectivity of the 
graph is low at some areas. For example, the lower right corner blob is only singly-connected to the rest of the graph. Some (blue) start locations may overlap 
with some (red) goal locations. 


computational results, each data point in the figures is an average over 10 sequentially, randomly generated instances. For each 
obstacle percentage, we start from the lowest number of robots (usually 10 or 20) and allocate a maximum of 600 seconds for 
each problem instance. If an instance takes more than 600 seconds to produce a result, the instance is failed and we move to 
the next obstacle percentage. This also means that a data point is given only if each of the 10 instances is completed within 
600 seconds. Over the same set of problem instances, the MinMakespan algorithm is executed in the exact manner (which 
produces optimal makespan solutions) and with the k-way split heuristic. 

The exact makespan computation result is summarized in Fig. 


10 For all obstacle settings, the MinMakespan algorithm 


computes optimal makespan solutions consistently for up to 100 robots with an average computation time of no more than 
100 seconds, notably, for 50 robots and obstacles up to 20%, the MinMakespan algorithm is able to complete in about 
15 seconds in all cases. From the top plot of Fig. 10 we observe that for each fixed obstacle percentage, the computation 


time appears to grow exponentially with respect to the number of robots. The computational difficulty of a particular problem 
instance also heavily depends on the actual optimal makespan. For example, a problem instance in the case of 25% and 20 
robots has a particularly long makespan (see the bottom plot of Fig. [T^, resulting an unexpected jump of the computation 
time. 

Whereas the exact MinMakespan algorithm is reasonably efficient, the k-way split heuristic brings a significant per¬ 
formance boost, allowing a much higher robot-vertex density in general. In our tests, we are often able to double or triple 
the supported robot-vertex density. For the 24 x 18 grid, we evaluated the k-way split heuristic for k up to 16. The 4-way 
split performance is illustrated in Fig. 11 In the figure, we measure optimality using a conservatively estimated optimality 


ratio. To obtain this, we divide the objective value returned by the optimizer over a conservative estimate (a lower bound). 
For makespan, this lower bound estimate is obtained by first computing the shortest path for each robot ignoring all other 
robots. The minimum makespan estimate is obtained by taking the maximum length over all these shortest paths. Clearly, the 
optimality ratio obtained this way is an underestimate. 

We make three comments over Fig. 11 First, from the top plot, we observe that our method is highly effective in terms of 
computation time, capable of computing minimum makespan solutions for up to 180 robots, which translates to a maximum 
robot-vertex ratio of 44%. The majority of the cases are solved in under 10 seconds. Even when there are 25% obstacles, 
we could solve the problem consistently for up to 60 robots in about 40 seconds. Second, we again observe an exponential 
relationship between computation time and the number of robots. Third, all computed solutions are very close to being optimal, 
with all but one case having an optimality ratio of below 1.1. The average minimum makespan for these instances is about 35. 
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Number of Robots 



Number of Robots 


Fig. 10. [top] Average computation time of the exact MinMakespan algorithm over instances on a 24 x 18 grid with randomly placed obstacles and 
start/goal locations, [bottom] The (average) optimal makespan. 




Number of Robots 


Fig. 11. [top] Average computation time of MinMakespan algorithm (with 4-way split heuristic) over instances on a 24 x 18 grid with randomly placed 
obstacles and start/goal locations, [bottom] The achieved (conservatively estimated) optimality ratio. 


The rest of the fc-way split evaluation is presented in Fig. 12 in which the computation time and optimality ratio are shown 
side by side without the axis labels. We also omit the key of the plots, which is the same as those from Fig. representing 
different obstacle percentages. In conjunction with Fig. 10 and Fig. 11 as fc increases, we observe a general trend of reduced 
computation time at the expense of some optimality loss. With 16-split, MinMakespan can solve problems with 300 robots, 
corresponding to a robot-vertex ratio of 69%. 

For comparison, we ran ID-based anytime algorithm over the same set of instances with a 600-second time limit (we also 
attempted ODh-ID and WHCA*, which do not go past 40 robots under the same setup). The result is plotted in Fig. 13 
ID actually performs quite well for up to 60 robots, which can be attributed to its A* root with minimum overhead as 
compared to our method. However, the performance of ID degrades faster-it does not scale well beyond 100 robots in our 
tests. MinMakespan, with 2-way split, readily outperforms ID when there are 40 or more robots. Conceivably, it may be 
possible to combine ID and fc-way split to make it run faster. However, adding fc-way split to ID will inevitably make the 
overall makespan more sub-optimal. 
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Fig. 12. Performance of the MinMakespan algorithm with fc-way split for k = 2 (top row), 8 (middle row), and 16 (bottom row). The computation 
time and optimality ratio plots for each k are shown side-by-side. The a:-axis for all plots represents the number of robots. The y-axis for the plots on the 
left represents the computation time in seconds. The y-axis for the plots on the right represents the optimality ratio. The keys for all plots are the obstacle 
percentage, identical to that of Fig. |10| 




Fig. 13. Performance of the ID-based anytime algorithm over the same set of problem instances. The plot setup is the same as that from Fig.|12| 


B. Minimum Makespan Solution to -puzzles 

Next, we evaluate the efficiency of the algorithm MinMakespan for finding minimum makespan solutions to the V^-puzzle 
for n = 3,4, 5, and 6. These problems have a robot-vertex ratio of 100%, making them highly constrained and extremely 
challenging. Note that an V^-puzzle instance is always solvable for n > 3 (see the Appendix); this means that all states are 
connected in the state (search) space. We ran Algorithm MinMakespan on 100 randomly generated TV^-puzzle instances for 
n = 3,4,5. For the 9-puzzle, computation on all instances completed successfully with an average computation time of 0.46 
seconds per instance. To compare the computational result, we implemented a (optimal) BFS algorithm. The BFS algorithm is 
heavily optimized: For example, cycles of the grid are precomputed and hard coded to save computation time. Since the state 
space of the 9-puzzle is small, the BFS algorithm is capable of optimally solving the same set of 9-puzzle instances with an 
average computation time of about 1.08 seconds per instance. 

Once we move to the 16-puzzle, the power of general ILP solvers becomes evident. MinMakespan solved all 100 randomly 
generated 16-puzzle instances with an average computation time of 4.2 seconds. On the other hand, the BFS algorithm with 
a priority queue that worked for the 9-puzzle ran out of memory after a few minutes. As our result shows that an optimal 
solution for the 16-puzzle generally requires 6 time steps, it seems natural to also try bidirectional search, which cuts down the 














14 


total number of states stored in memory. To complete such a search, one side of the bidirectional search generally must reach 
a depth of 3, which requires storing over 5 x 10® states (the branching factor is over 1000), each taking 64 bits of memory. 
This translates into over 4GB of raw memory and over 8GB of working memory, which is more than the JavaVM can handle: 
A bidirectional search ran out of memory after about 10 minutes in general. We also experimented with C++ implementation 
using STL libraries, which yields similar result {i.e., ran out of memory before reaching a search depth of 3). 



Fig. 14. An instance of a 25-puzzle problem solved by MinMakespan. 


For the 25-puzzle, without a good heuristic, bidirectional search cannot explore a tiny fraction of the fully connected state 
space with about 10^® states. On the other hand, MinMakespan again consistently solves the 25-puzzle, with an average 
computational time of 391.6 seconds over 100 randomly created problems. Fig. 14 shows one of the solved instances with a 
7-step solution given in Fig. 15 Note that 7 steps is obviously the least possible since it takes at least 7 steps to move robot 10 
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Fig. 15. An optimal 7-step solution (from left to right, then top to bottom) to the 25-puzzle problem from Fig.[T^ by MINMAKESPAN in about 15 minutes. 

to its desired goal. We also briefly tested MinMakespan on the 36-puzzle. While we had some success here, MinMakespan 
generally does not seem to solve a randomly generated instance of the 36-puzzle within 24 hours, which has 3.7 X 10^^ states 
and a branching factor of well over 10®. 

As comparison, WHCA* can not solve 9-puzzle. OD+ID, and ID and can only solve 9-puzzle consistently and cannot 
solve any 16-puzzles in 600 seconds. 


C. Minimum Makespan on 8 x 8, 16 x 16, and 32 x 32 Grids 

In this section, we evaluate MinMakespan with underlying graphs that are 8x8 grids, 16 x 16 grids, and 32 x 32 grids 
with 20% obstacles. In addition to further demonstrating the effectiveness of MinMakespan, this allows us to better compare 
our results. 8x8 and 16 x 16 grids are used as the environment for evaluation in Surynek (20121. 32 x 32 grids with 20% 
obstacles are used for evaluation in [Silver (2005|l; Standley and Korf (201 l|l. 

For 8x8 and 16 x 16 grids, the instances are constructed using the same procedure stated in Section VI-A| Again, each 
data point is an average over 10 sequentially randomly created instances. Given the size of 8 x 8 and 16 x 16 grids, we 
limit k to 8; using 16-way split can solve more instances but incur an average solution optimality between 2-4 optimal. Each 
instance is given a time limit of 600 seconds. The outcome of these experiments is plotted in Fig. 16 along with the result 
from running ID. OD+ID and WHCA* cannot consistently solve the instances with 20 robots in 600 seconds. Makespan, 
instead of optimality ratio, is used in Fig. 16 for easy comparison with the results from [Surynek ( 2012| l. 

We observe that, over the 8x8 grid, with 2-way split heuristic, MinMakespan can solve problems with 50 robots to 
almost true optimal solutions in just 10 seconds. With 4-way split, MinMakespan can further push to 60 robots (robot-vertex 
ratio of 94%) with solutions that are within 1.7-optimal. ID can only handle up to 30 robots. As reported in Surynek (2012 1 , 
COBOPT generally takes more than half an hour to produce its final solution when there are 24 or more robots]^ The solution 

^We did not directly run COBOPT over our randomly created instances. However, the instances in | Surynek] {2012) are created in an identical manner. 
Therefore, given that similarly powered computers are used, the computation time and solution makespan are directly comparable between ours and those 
from [Surynek] {2012). 
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Fig. 16. Performance of the MinMakespan algorithm with k-way split for k = 2-8 and the ID-based anytime algorithm. The computation time and 
solution makespan plots are shown side-by-side. The x-axis for all plots represents the number of robots. The y-axis for the plots on the left represents the 
computation time in seconds. The y-axis for the plots on the right represents the solution makespan. Note that we did not plot the makespan computed by the 
exact MinMakespan algorithm. Instead, the lower bound (LB) estimate of makespan is plotted (in magenta color), [top] Result on the 8x8 grid, [bottom] 
Result on the 16 x 16 grid. 


quality also degrades quickly as the number of robots increases. For example (Fig. 2 and Fig. 3 in [Surynek ( |2012| l), at 50 
robots, COBOPT takes over an hour to produce a solution with a makespan of over 160, whereas our 2-way split heuristic 
yields a near-optimal makespan of 11.8 in just 10 seconds. 

Over the 16 x 16 grid, MinMakespan is able to handle instances with 190 (robot-vertex ratio of 74%) robots with 8-way 
split while at the same yielding solutions that are always less than 1.4-optimal. When switched to 4-way split, MinMakespan 
can consistently solve problems with up to 160 robots to no more than 1.03-optimal. In comparison, ID can solve instances 
with up to 80 robots to relatively good quality. Taking on average an hour of computation, COBOPT can handle up to 128 
robots; the solution quality is quite poor. For example (Fig. 4 and Fig. 5 in |Suryn^ ( |2012| l), for about 100 robots, the computed 
makespan by COBOPT is at 200 whereas the optimal makespan is about 25. Across 8x8 and 16 x 16 grids, we generally 
observe a speedup of over 100 times when MinMakespan (with the 4-way split heuristic) is compared with Cobopt. At 
the same time, our method yields solutions with much smaller makespan. 

A classical test scenario is the 4-connected 32 x 32 grid with 20% vertices randomly removedfjFor completeness, we also 
perform a brief evaluation of this setup. We randomly generate the instance as before, run the test, and plot the result in 
Fig. 17 From the figure, we observe a pattern consistent with experiments on the 8 x 8, 16 x 16 and 24 x 18 grids. 




Fig. 17. Performance of MINMAKESPAN and the ID-based anytime algorithm over 32 X 32 grid. The axis setup is the same as that from Fig. |12| 


D. Algorithm Performance over All Objectives 

Last, we evaluate the performance of our algorithms at optimizing all Objectives [T]|^ The result on MinMakespan is already 
presented in Section VI-A| which we also use as the solution to optimizing Objective (/.e., we simply use MinMakespan 


®Some work (e.g., Standley and Korf] 2011 ) also adopts an 8-connected model. That is, each vertex is on the grid is connected to its 8-neighborhood. 
This causes unit cost to be assigned to all edges, although a diagonal edge should have length \/2 times that of a non-diagonal edge. Since we are modeling 
robots in this work, we do not discuss the 8-connected model here. However, we mention that our algorithms easily extend to 8-connected model. Our tests 
show that we can in fact compute near-optimal makespan for 400 robots on 32 x 32 grids assuming 8-connectivity. 
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in place of MinMaxDist due to its superior performance). Note that no change to the plots are needed here because, on 
one hand, we can use a (near)-optimal makespan solution as a solution to minimizing maximum single-robot traveled distance. 
This is true because the makespan of a solution is always no less than the minimum maximum single-robot distance. On the 
other hand, the lower bound estimate for minimum makespan is the same as that for minimum maximum single-robot distance. 

Before moving to MinTotalTime and MinTotalDist, we note that these algorithms possess some properties of an 
anytime algorithm, which is of practical importance. In solving these ILP models, the solver generally uses variations of the 
branch-and-bound algorithm Land and Doig ( 1960| l. For computing total time (and distance) optimal solutions, a branch-and- 
bound algorithm always try to find a feasible solution first and then iteratively improve over the feasible solution. This naturally 
leads to improving solution quality commonly observed in an anytime algorithm. The anytime property of MinTotalTime 
and MinTotalDist allows us to set a desired sub-optimal threshold to reduce the computation time. Note that the same 
cannot be said for MinMakespan because a feasible solution here is the optimal solution. 

Our next set of results focuses on the MinTotalTime algorithm (Fig. [TSl l. The general setup is the same as that used in 
Section VI-A In particular, the same set of problem instances is used. In out experiment, we limit both time (600 seconds) 
and required sub-optimality threshold (automatically adjusted) to achieve a balanced performance. The lower bound estimate 
for computing optimality ratio is obtained by summing over the individual shortest path lengths. The actual optimal total time 
is about 15 times the number of robots {i.e., 150 for 10 robots and 1500 for 100 robots), regardless of the percentage of 
obstacles. We observe that the MinTotalTime algorithm is fairly effective, capable of computing 1.1-optimal solutions for 
up to 100 robots in the allocated time. 

Similar outcome is also observed in the performance evaluation of the MinTotalDist algorithm with the 4-way split 
heuristic (Fig. [T9| ). The optimal total distance is again about 15 times the number of robots. In comparison with the total time 
optimal case, due to the 4-way split heuristic, the MinTotalDist algorithm is faster but produced solutions that are more 
sub-optimal but still quite good. 
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Fig. 18. [top] Average computation time of MinTotalTime algorithm over instances on a 24 x 18 grid with randomly placed obstacles and start/goal 
locations, [bottom] The achieved (conservatively estimated) optimality ratio. 


For comparison, we run MinTotalTime, MinTotalDist, and ID (total time and total distance versions) on 24 x 18 
grids with 0-20% obstacles with a maximum time limit of 600 seconds. The setting is slightly different from that used in 
obtaining Fig. [T^ and 19 we do not set a sub-optimal threshold here. The result is plotted in Fig. 20 In the case of total time 
(Objective [^, ID could solve more instances. For the instances that are solved, the achieved optimality is similar. For total 
distance (Objective]^, we observe similar outcome. Here, ID could produce one more data point; the achieved optimality by 
both methods are again comparable (and good). Overall, MinTotalTime, and MinTotalDist are competitive with ID. 


VII. Conclusion 

In this paper, we propose a general algorithmic framework for solving MPP problems optimally or near-optimally. Through an 
equivalence relationship between MPP and network flow, we provide ILP model-based algorithms for minimizing the makespan 
(last arrival time), the maximum (single-robot traveled) distance, the total arrival time, and the total distance. In conjunction with 
additional heuristics, our algorithmic solutions are highly effective, capable of computing near-optimal solutions for hundreds 
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Fig. 19. [top] Average computation time of MinTotalDist algorithm over instances on a 24 x 18 grid with randomly placed obstacles and start/goal 
locations, [bottom] The achieved (conservatively estimated) optimality ratio. 






Fig. 20. Performance of MinTotalTime MinTotalDist, and the ID-based anytime algorithm over 24 x 18 with 0-20% obstacles. Each instance is 
allowed 600 seconds of time. The axis setup is the same as that from Fig. |12| [top] MinTotalTime versus total time ID; the red lines correspond to data 
for MinTotalTime. [bottom] MinTotalDist versus total distance ID; the red lines correspond to data for MinTotalDist. 


of robots in seconds in scenarios with very high robot-vertex ratio. In pushing for high performance algorithms aim at solving 
MPP optimally or near-optimally, our eventual goal is to apply these algorithms to multi-robot path planning problems in 
continuous domains-our preliminary work toward this goal has begun to show promising results, producing algorithms that 
can compute solutions for around a hundred disc robots in 2D environments with holes. 
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Appendix 

A. Properties of the N^-puzzle 

The example problem from Fig. [I] easily extends to an V x V grid; we call this class of problems the -puzzle. Such 
problems are highly coupled: No robot can move without at least three other robots moving at the same time. At each step, all 
robots that move must move synchronously in the same direction (per cycle) on one or more disjoint cycles (see e.g.. Fig. [^. 
To put into perspective the computational results on V^-puzzles that follow, we make a characterization of the state structure 
of the V^-puzzle for N > 3 (the case of V = 2 is trivial). 



Fig. 21. A 3-step procedure for exchanging robots 8 and 9. 


Proposition 6 All states of the 9-puzzle are connected via legal moves. 


Proof. We show that any state of a 9-puzzle can be moved into the state shown in Fig. [TJb). From any state, robot 5 can be 
easily moved into the center of the grid. We are left to show that we can exchange two robots on the border without affecting 
other robots. This is possible due to the procedure illustrated in Fig. 21 □ 


Larger puzzles can be solved recursively: We may first solve the top and right side of the puzzle and then the left over 


smaller square puzzle. For a 16-puzzle, Fig. 22 outlines the procedure, consisting of six main steps: 

1) Move robots 1 and 2 to their respective goal locations, one robot at a time (first 1, then 2). 

2) Move robots 3 and 4 (first 3, then 4) to the lower left corner (top-middle figure in Fig. [2^ . 

3) Move robots 3 and 4 to their goal location together via counterclockwise rotation along the cycle indicated in the 


top-middle figure in Fig. 22 


4) Move robot 8 to its goal location. 

5) Move robots 12 and then 16 to the lower left comer. 

6) Rotate robots 12 and 16 to their goal locations. 

It is straightforward to see that larger puzzles can be solved similarly. We have thus outlined the essential steps for proving 
Proposition l^below; a more generic proof can be written using generators of permutation groups, which we omit here due to its 
length. Proposition implies that, for N > 3, all instances of -puzzles are solvable. The constructive proofs of Proposition 
1^ and 1^ lead to recursive algorithms for solving any A^^-puzzle (clearly, the solution is not time/distance optimal in general). 


Theorem 7 All states of an N^-puzzle, N >3 are connected via legal moves. 


Corollary 8 All instances of the -puzzle, N > 3, are solvable. 

By Proposition]^ since all states of a -puzzle for A^ > 3 are connected via legal moves, the state space of searching an 
W^-puzzle equals factorial. For 16-puzzle and 25-puzzle, 16! > 10^^,251 > 10^®. Large state space is one of the three 
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Fig. 22. A solution scheme for solving top/left sides of the 16-puzzle. 


reasons that make finding a time optimal solution to the -puzzle a difficult problem. The second difficulty comes from the 
large branching factor at each step. For a 9-puzzle, there are 13 unique cycles, yielding a branching factor of 26 (clockwise and 
counterclockwise rotations). For the 16-puzzle, the branching factor is around 500. This number balloons to over 10^ for the 
25-puzzle. This suggests that on typical commodity personal computer hardware (assuming a IGHz processor), a baisc breadth 
first search algorithm will not be able to go beyond depth of 3 for the 16-puzzle and depth 2 for the 25-puzzle in reasonable 
amount of time. Moreover, enumerating these cycles is a non-trivial task. The third difficulty is the lack of obvious heuristics: 
Manhattan distances of robots to their respective goals prove to be a bad one. For example, given the initial configuration as 
that in Fig. [Ja), the first step in the optimal plan from Fig. |^gets robots 1, 3, 4, 6, 8, 9 closer to their respective goals while 
moving robots 2, 7 farther. On the other hand, rotating counterclockwise along the outer cycle takes robots 1, 3, 4, 5, 6, 8, 9 
closer and only moves robot 7 farther. However, if we instead take this latter first step, the optimal plan afterward will take 5 
more steps. 


















































